The Dynamical Implications of Multiple Stellar Formation Events 

in Galactic Globular Clusters 



Jonathan M.B. Downing 1 and Alison Sills 

Department of Physics and Astronomy, ABB-2^1, McMaster University, 1280 Main Street 

West, Hamilton, Ontario, L8S J^Ml, Canada 

downin@ari .uni-heidelberg. de , asills@mcmaster . ca 

<N ■ 

X) '■ ABSTRACT 

<D ■ 

plh : 

^ \ Various galactic globular clusters display abundance anomalies that affect the 

CN \ morphology of their colour-magnitude diagrams. In this paper we consider the 

\ possibility of helium enhancement in the anomalous horizontal branch of NGC 

2808. We examine the dynamics of a self-enrichment scenario in which an initial 
7-h ! generation of stars with a top-heavy initial mass function enriches the interstellar 

^ \ medium with helium via the low-velocity ejecta of its asymptotic giant branch 

O stars. This enriched medium then produces a second generation of stars which 

are themselves helium-enriched. We use a direct N-body approach to perform 
five simulations and conclude that such two-generation clusters are both possible 
P-f and would not differ significantly from their single-generation counterparts on 

the basis of dynamics. We find, however, that the stellar populations of such 
clusters would differ from single-generation clusters with a standard initial mass 
function and in particular would be enhanced in white dwarf stars. We conclude, 
. ^ | at least from the standpoint of dynamics, that two-generation globular clusters 

are feasible. 

Subject headings: stellar dynamics — methods: N-body simulations — globu- 
lar clusters: general — globular clusters: self-enrichment — globular clusters: 
individual(NGC 2808) — stars: AGB — stars: HB 



1. Introduction 

Globular clusters are gravitationally bound collections of stars, typically including be- 
tween 10 3 and 10 7 members, which normally occur within the halos of galaxies. Galactic 
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globu lar clusters are understood to be old objects, up to 12 Gyr in age (IMeissner h Weiss 
20061 ). and are thus primordial components of the galaxy. In the simple picture, all stars 
within a globular cluster are thought to have formed at the same time and out of the same 
medium, thus sharing the same age and chemical composition. In particular all stars in a 
given globular cluster have the same value of [Fe/H] and are said to be mono-metallic. The 
only differentiation between stars in this scenario is their spectrum of masses determined 
by an initial mass function (IMF). There is thought to be a universa l IMF for all galactic 
globular clusters but there is some debate as to its exact form (jKroupal (120021 ) and references 
therein) . 

Although this description works well for the most part, several globular clusters display 
population characteristics that are inconsistent with this model. In many clusters all stars 
have the same [Fe/H] but there are star-to-star abundances variations in other elements. In 
partic ular many globular c lusters stars show strong star-to star O-Na and Mg-Al anticorrela- 
tions (IGratton et al.ll2004j ). Some clusters also display peculiar morphologies in the horizon- 
tal branch of the Colour-Magnitude Diagram (CMD). In some cases different clusters with 
the same [Fe/H] and with otherwise identical colour-magnitu de diagrams will have differing 
horizo n tal branch morpho logies (e.g. the case of M3 and M13 (ICatelan fc de Freitas Pacheco 
(119951 ) . iKraft et al.l (Il992l )). This is called the "second parameter problem" since such clus- 
ters seem to require a second parameter (other than metalicity) in order to explain the vari- 
ations between them. In some cases, the presence of very extended horizontal branch blue 
tails is linked to dynamical or abundance effects. One possible source of these anomalies is 
a self-enrichment scenario where chemicals produced by some cluster stars are incorporated 
into other cluster stars which then display unusual characteristics such as enhancements 
of light and heavy elements (particularly products of hot-hydrogen burning and the CNO 
cycle) and unusual c olours (frequently bluer) with respect to field stars (see the review by 
Gratton et al.l (12004 ) for further details). We propose to explore the dynamics of a particular 
variation of the self-enrichment scenario in which excess helium is produced by asymptotic 
giant branch (AGB) stars and which may generate the anomalous horizontal branch of the 
globular cluster NGC 2808. 

NGC 280 8 is a southern galactic globular cluster first observed in detail in the late 1960s 
(|ALcainolll96sh . . t is quite massive at 1. 6 x 1O 6 M , has a moderately high velocity dispersion 
a n = 13 Akm /s (jPryor fc Meylanl 119931 ) and is relatively metal rich with [Fe/H] = —1.09 
dHarrislll996h . As early as 1974 a peculiar morphology (both bimodal and extende d) was 
discovered in the horizontal branch (HB) of the CMD of NGC 2808 Jflarrisl [19741 ) . The 
horizontal branch of NGC 2808 consists of both a small horizontal clump (hearafter called 
the red horizontal branch or RHB) and a separate, extended tail reaching vertically from 
16th to 22nd magnitude (hereafter called the blue horizontal branch or BHB). The BHB 
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is also multi-modal with three distinct groups (called from top to bottom extended blue 
tails (EBTs) 1, 2 and 3). This is very different from canonical HB morphology (flat and 
monomodal) and in order to explain it a special scenario is required. 

One way to explain the anomalous morphology of the HB of NGC 2808 is to assume 
that the BHB stars are he lium enriched while th e RHB stars contain the cosmological helium 
abundance of Y ~ 0.24 jD'Antona et alJl2002h . The enrichment cannot be attributed to 
helium variations across the primordial cloud s ince the necessary levels of enrichment are 
extre me (Y ~ 0.32 - 0.4 is needed in EBT 3 JD'Antona fc CaloJliooi ). JD'Antona et al. 



20051 )) and the helium must arise from stellar processes through self-enrichment. It has been 
argued (IGnedin et al.ll2002l ) that the more massive globular clusters (particularly cu-Centari) 
may have sufficient self-gravity to retain the low-velocity ejecta from AGB stars (and in the 
case of the most massive even some of the ejecta from Type II supernova). In the case of 
clusters within ~10 kpc (such as cu-Centari) IGnedin et al.l (120021 ) argue that interaction with 
the galactic disc would strip such ejecta from clusters befor e it can expe rience star formation. 
NGC 2808, however, is 11.1 kpc from the galactic center (IHarrisI 1 1 9961 ) a nd may avoid such 
stripp ing long enough to undergo a second generation of star formation (jD'Antona &: Caloi 
20041 ) . This leads to a scenario in which helium self-enrichment from AGB stars produces 
the horizontal branch morphology of NGC 2808. In particular iD'Antona fc Caloil (120041 ) 
propose a model where NGC 2808 has experienced two stellar formation events. The first 
event features a top-heavy mass function that enhances the number of 3 — 5M stars. These 
evolve to the asymptotic giant branch (AGB) phase in < 200 Myr and enrich the interstellar 
medium (ISM) with helium (and other elements) through stellar winds. A second generation 
of stars then forms from the helium-enriched ISM. An age difference of < 200 Myr would 
not be observable because it is a small fraction of the lifetime of a low mass star, but the 
helium difference in the populations would be. The low-mass stars in the first generation 
would form the RHB of NGC 2808 and the enriched second generation would form the 
BHB. The multi-modality of the BHB (EBTs 1, 2 and 3) can be explained by varying levels 
of helium enrichment produced either by several minor stellar formation events producing 
different levels of helium enrichment or possibly by differential helium enrichment across the 
cluster. Modeling more than two generations or trying to account for spatially dependent 
helium enhancement is beyond the scope of our simulations and in this paper we do not 
investigate the substructure in the BHB. It is possible that helium enrichment is the source 
of other anomalous CMD morphologies. Notably the double main sequence of the massive 
cluster uj Centauri can be e xplained if the b l ue m ain sequence stars are helium enriched 
( Piotto et al. 2005 1"*! Indeed D'Antona et al. ( 2005 ) claim helium enrichment in some 20% 



1 w Centauri has many other peculiarities and may be the core of a stripped dwarf galaxy. It cannot be 
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of main sequence stars for NGC 2808 (although this case is less clear cut) but we have not 
tried to model this effect. 



Based on stellar evolution models iD'Antona fc Caloil (120041 ) claim that, given certain 
assumptions, sufficient helium can be produced by AGB stars in the two-generation scenario 
to explain the BHB of NGC 2808. If enhancements in other elements produced by AGB 
stars are ignored and the top-heavy mass function is accepted, the scenario may be plausible 
on the basis of stellar chemistry. The massive first generation may, however, have effects on 
both the dynamics and population of the cluster. Specifically the high mass of the initial 
generation and the effect of adding a second generation to a dynamically evolved object 
could lead to dynamical instabilities or affect the final spatial distribution of a star cluster. 
In addition the first generation will leave many intermediate- mass remnants , prima rily white 
dwarfs, which should still be observable after 10-12 Gyr (e.g. lHansen et al.l (120021 ) ). We use 
a direct N-body approach to explore the dynamical evolution of a globular cluster in the 
D'Antona and Caloi two-generation formation scenario. 



2. Method 



To perform our simulations we use Starlab, a dire ct N-body evolu tion code. Starlab 
uses a fourth-order Hermite predictor-corrector scheme ( McMillan 1986h and contains a fu ll 
suite of stellar and binary evolution al gorithms based on the work o f lEggelton et al.l (119891 ). 



Starlab is described in some detail in 



Portegies Zwart et al.l (120011 ) and is freely available 



with some documentation from http: //www.ids.ias.edu/~sta rlab/. We also use a GRAPE- 
6A N-body gravitational accelerator (IFukushige et al.ll2005l ) to improve the speed and size 
of our simulations. 



2.1. Top-Heavy IMF 



The D'Antona and Caloi top-heavy IMF given in lD'Antona fc Caloil (120041 ) is a broken 
power-law which takes the form: 



dN 
dM 



Cl M-( 1+Q ) if M < M B 
c 2 M^ 1+ ^ if M > M B 



and is defined by the exponents a and j3 and a transition mass Mb {c\ and c 2 are normaliza- 
tion constants). The explicit implementation of this IMF in Starlab is discussed in Appendix 



considered typical of galactic globular clusters 
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A. A D'Antona and Caloi IMF created with this implementation where a = —0.5, (3 = 3.0 
and Mb = 3.8 is compared with a standard Sa l peter IMF in Figure HJ These parameters 
correspond to the IMF from lD'Antona fc Caloil (120041 ) with the lowest overall mass still ca- 
pable of providing the requisite helium enhancement. This is ostensibly the least extreme 
scenario and, since discrepancies between simulations of single-generation clusters and more 
extreme examples of the D'Antona and Caloi scenario might not preclude agreement in less 
extreme cases, we use these parameters in all of our simulations. 



2.2. Details of the First Generation 

In modeling the first generation there are three compontents to consider, one stellar and 
two gaseous. The stellar component is simply the stars of the first generation and can be 
dealt with by standard N-body modeling. The first gas component consists of the remnant 
gas from the first generation of star formation. The second gas component is composed of the 
ejecta from the AGB stars of the first generation which produces the helium enhancement 
for the second generation. 

There is little information available on the configuration of remnant gas (first compo- 
nent) in young globular clusters and for simplicity we model the remnant ISM left over from 
the first generation of star formation as an analytic Plummer potential, 0(r) = — y=§=j=g ■ 
This potential is not static since the remnant ISM is stripped from the young cluster by the 
radiation from young O stars, early supernovae and possibly due to interactions with the 
galactic disk. The exa ct rate of g as str i pping from young clu s ters is unknown but it on the 



order of 10s of Myrs (jGoodwinl (119971 ). iBastian fc Goodwin! (120061 )). It is clear, however 



that this mass-loss affects early cluster evolution since the potential in which the cluster 
lives becomes more shallow over time and the cluster can expand and in the extreme case 
be completely disrupted. 

The stellar ejecta forming the second compone nt originates p r imari ly from supernovae 



and the stellar winds of AGB stars. According to iGnedin et al.l (120021 ) the characteristic 
terminal velocity for AGB winds is ~ 15 km s -1 whereas the escape velocity for a cluster of 
the mass of NGC 2808 is on the order of 60 km s -1 in the core and 40 km s _1 at the half 
mass radius. It is plausible that a large fraction of the mass lost by AGB stars will in fact be 
retained as gas by the cluster. In our models we assume that all of the mass lost due to stellar 
evolution is retained, that is, we do not treat ejecta from different progenitor stars separetly. 
This assumption is not necessarily realistic since the high-velocity supernovae ejecta may 
be able to escape from the potential well of the young cluster and some of the AGB ejecta 
may be stripped by radiation in the same way as the remnant ISM. The majority of stars 
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in the first generation are below 1OM and will not undergo supernovae events in the first 
200 Myr of the cluster's life. Therefore the AGB ejecta should be the dominant component 
and including the supernovae ejecta is not a concern. The stripping fraction for the AGB 
ejecta is also not clear and treating it would vastly expand our parameter space without 
providing any particular physical insight. Thus the choice of a 100% re tention fraction for 



the el ecta is a reasonable first approximation and is also in the spirit of iD'Antona fc Caloi 
(120041 ) . We model the retained ejecta as another single, global Plummer potential with a 
mass increasing at the rate mass is lost from the stars due to stellar evolution. To be more 
realistic, we should add gas to the simulation locally where it is lost by the stars. Since at 
\bkms~ 1 the gas crossing time at the initial viral radius is ~ 3 Myr (about two simulation 
time units) and the particle crossing time is ~ 8.5 Myr (about four simulation time units) 
we expect gas produced locally to be quickly distributed throughout the cluster. Thus the 
assumption of a global profile for the second gas component seems reasonable. 

Since Starlab does not contain prescriptions for time-evolving potentials and to incor- 
porate such prescriptions linked both to radiation pressure and stellar evolution would be 
time-consuming we make the following approximations: (1) We assume that after 200 Myr 
all of the remnant ISM has been removed from the cluster. (2) We assume that the remnant 
ISM is stripped from the cluster at the same rate at which ejecta from the first generation 
stars is added to the cluster. (3) We assume that the initial mass of the first component 
is the same as the final mass of the second component. This leads to a model where the 
first-generation is embedded in a static Plummer potential with a total mass equal to the 
mass lost by stellar evolution in the first generation. Assumption (1) is almost certainly 
secure since young clusters appear to be stripped of their remnant star forming material 
on timescales of a few 10s of Megayears but this means assumption (2) is not necessarily 
accurate since the remnant ISM will be ejected more rapidly than it is replaced by the ejecta 
from the first generation (added to the cluster over a timescale of 200 Myr). We argue that 
the two interacting potentials will act in three phases. Initially the remnant ISM potential 
will provide a deep well which will cause the cluster to contract. Once the remnant ISM 
has been removed and before the ejecta potential becomes large the cluster will be able to 
expand normally due to mass-loss. Finally the ejecta potential will dominate and strongly 
mitigate the effect of mass-loss and again slow the expansion. The global static potential 
will have the correct behaviour at the beginning and end and will simply slightly deepen 
the potential in the middle phase. In both cases the dynamical effect should simply be to 
suppress the expansion of the cluster due to stellar-evolution mass-loss. Assumption (3) is 
simply convenient and, since the exact star formation rate in young clusters is unknown, 
any variation of this paramater would not provide any new physical insight. To confirm the 
static potential has no additional dynamical effects, we compare the evolution of one of our 
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first generations with and without a potential in Figure [2j 

As expected the only qualitative difference between the two plots in Figure [2] is that 
with the addition of the Plummer potential both the massive and overall Lagrangian radii 
increase more slowly and lead to a smaller, more compact final cluster. This effect is seen in 
both the massive and overall Lagrangian radii and is a good indication that mass segregation 
has not been suppressed by the addition of the exte rnal field. Since mass-segregation is a 



two-body relaxation effect (IFarouki fc Salpeterl Il982l ) this shows that two-body relaxation 



also has not been affected by the addition of the Plummer potential. 

Obviously a pair of time-evolving potentials linked to the stellar physics of the code 
would be preferable but the difference would be in details rather than in overall effect. 
Any galactic globular cluster will have experienced ~ 10 Gyr of dynamical evolution since 
the potential formation of a second generation and thus be a fully dynamically relaxed 
object. Since dynamically relaxed N-body systems are insensitive to their initial conditions 



( iBinney fc Tremaind 119871 ) we would not expect fine details of the gas potential in the first 
generation to be significant for current observations. We note that our predictions are not 
applicable to two-generation clusters in the process of forming their second generation or 
even those less than a relaxation time old. For such objects, however, it is not clear that 
global potentials would be sufficient to accurately resolve the details of evolution and such 
models may have to wait until hybrid SPH-N-body codes are available. 



2.3. The Second Generation 



In order to match current observations we assume the secon d generation fo rms with a 
more standard IMF and for simplicity we choose a Salpeter IMF (jSalpeterl Il955h . In order 
to implement the second generation we must find some way add stars to the output file 
from the simulation of the first generation (in Starlab such a file is know as a snapshot). 
The easiest way to implement this in Starlab is to create a separate snapshot using the 
parameters characterizing the second generation and then combining this snapshot with 
the output snapshot from the first generation using the routine merge_snaps. We choose 
the total mass of the second-generation to match the mass of the Plummer potential and 
hence the mass lost in our first-generation. This implicitly im plies 100% star form ation 



impi 

efficiency. This assumption is almost certainly unrealistic (e.g. lLada fc Ladal (120031)) but 



it is both in keeping with the spirit of the scenario proposed by iD'Antona &: Caloil (120041 ) 
and it allows us to conserve mass when adding the second generation (all of the Plummer 
potential mass is replaced by stellar mass). Since the actual star formation efficiency during 
galactic globular cluster formation is not known, considering different star formation rates 
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would again increase our parameter space without providing any particular physical insight. 
We use the merge_snaps routine to combine the first and second generation snapshots into a 
single snapshot and manually remove the Plummer potential from the first generation. This 
procedure represents the conversion of a continuous gas potential into a generation of stars 
with the same total mass in one timestep (~ 1.5 Myr). Once the second generation has 
been added and the Plummer potential removed, the combined cluster is evolved for ~ 12 
Gyr which corresponds to 15-16 initial half-mass relaxation times (Tr#) (jSpitzerl Il987t ) . In 
order to produce simulations lasting more than 2-4 initial half-mass relaxation times it is 
necessary to introduce a stripping radius. This is because newly-formed neutron stars are 
ejected from the cluster with a high velocity due to their initial kick. Once far enough away, 
the code seems to treat pairs of such neutron stars as binaries and the three-body interaction 
between these and the rest of the cluster eventually reduces the timestep to zero and halts 
the simulation. Removing stars at a large radius from the cluster eliminates this problem. 



3. Initial Conditions 

We perform five simulations: a single generation model with a Salpeter IMF (SP), a 
single generation model with a D'Antona and Caloi top-heavy IMF (DC), and three two- 
generation models with a D'Antona and Caloi IMF for the first gener ation and a S a lpeter 



IMF for the second generation (Ca, Cb, and Cc). In all cases we follow iHurley et al.l (120051 ) 
and use Plummer models for our initial conditions. In the two-generation models we allow 
the first generation to evolve for ~ 200 Myr before adding the second generation according 
to a Salpeter IMF. In all cases the initial virial radius is the length unit and is set to 5.0 
pc. The mass unit for each simulation is the total mass of the simulated cluster. The energy 
and time units are defined by starting all simulations in virial equilibrium. Parameters for 
these simulations are given in Table [lj The number of particles used in these simulations 
is a c omprimise between achiving an accurate number of particles and our hardware limita- 
tions. iFukushige fc Heggid (119951 ) have explored the behavior of clusters with different mass 
funcitons and find that there is little difference between simulations carried out with 8192 
and 16384 particles. Practically, we find that ~ 12000 particles is the maximum number 
for which we can achive a reasonable wall-clock time and since we are in the middle of the 
8000-16000 particle range, we do not expect our results to scale to strongly with particle 
number. All clusters are simulated in isolation in order to keep our parameter space small 
and to eliminate effects that are not due to the addition of the second generation. We must, 
however, introduce a stripping radius for the reasons discussed above. We find that by strip- 
ping at 100 initial viral radii we eliminate the problematic fast-moving neutron stars while 
removing few other stars from the simulations. 
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4. Results 

We consider results for the dynamical and final stellar populations of our simulations 
separetly. 

4.1. Dynamics 

We present the initial half-mass relaxation times and some of the final physical data 
for our systems in Table [2j Each cluster has lost a significant portion of its initial mass. 
Few particles are removed from the simulations and this mass loss is due to stellar evolution 
rather than loss of stars from the cluster. This mass-loss is particularly pronounced in DC 
with its extreme initial enhancement of high-mass stars and is also a major factor in the two- 
generation clusters with their top-heavy first-generation IMFs. We find the initial half-mass 
relaxation times of the two-generation clusters to be somewhat smaller than that of the single 
generation Salpeter IMF while the half-mass relaxation time of the single-generation cluster 
with the D'Antona and Caloi IMF is much shorter than any of the other cases. Since each 
cluster has the same initial virial radius and thus the same initial physical concentration, the 
differences are due to the initial total mass increasing from SP through the two-generation 
clusters to DC. It is worth noting that the combined clusters undergo a sufficient number of 
relaxation times to allow the possibility of core collapse by the end of our simulations. 

In Figure [3] we present a comparison of the Lagrangian radii for the runs SP, DC and 
Ca (Cb and Cc are similar to Ca). The evolution of the Lagrangian radii for SP is slightly 
different from Ca and DC. In DC and Ca the Lagrangian radii grow quickly at first and then 
flatten significantly, with the 75% radius still growing slowly and the 50% and 10% flattening 
(and even contracting in the case of the 10% and 50% massive radii of Ca). By contrast the 
early growth period is much less evident in SP and all radii continue to grow slowly (with 
the 10% radius flattening at the end). The difference is due to the high-mass stars in DC 
and the first generation of Ca. These high-mass stars evolve quickly and consequently lose a 
great deal of mass over a short timescale (the first generation of Ca loses ~ 30% of its mass 
in 200 Myr). This mass-loss leads to a rapid reduction in cluster potential, causing rapid 
cluster expansion and a consequent increase in the Lagrangian radii. Once the high-mass 
stars evolve away, mass-loss becomes dominated by low-mass stars which evolve much more 
slowly leading to a much reduced mass-loss rate and a subsequent flattening of the Lagrangian 
radii. By contrast SP is always dominated by low-mass stars and experiences slow, constant 
mass-loss throughout its life, consequently displaying slow, constant expansion. 

We compare the Lagrangian radii of the most massive 10% of stars (cutoff masses given 
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in Table [2]) to the Lagrangian radii for the entire system in Figure [3j We note that there 
is a wide separation between the massive and overall Lagrangian radii in SP and Ca while 
for much of the life of DC the massive and overall Lagrangain radii are similar. That 
the massive Lagrangian radii are smaller indicates the high-mass stars are more centerally 
concentrated than the low-mass ones. This is a manifestation of mass segregation which is 
a natural consequence of equipartition of e nergy in phase space and is expected in all N- 
body systems with a range in initial masses (IBinney fc Tremaindll987l ). The smaller spread 
between the overall and massive population in DC seems to indicate that mass segregation 
is less advanced in this simulation than in SP and Ca. This is somewhat unexpected since 
for an N-body s ystem the timescale for m ass segregation is on the order of the two-body 
relaxation time (IBinney fc Tremaind 119871 ). Since DC has the shortest half-mass relaxation 
time of our simulations we would expect mass segregation would be the most advanced in 
this model. We postulate that the larger spread and higher overall mass in the IMF of DC 
(Figure [I]) both renders the overall statistics of this simulation less sensitive the the effects 
of a few very high-mass stars in the central regions and produces a more smooth transition 
between the low-mass halo and the high-mass core. Both of these processes could bring the 
overall and massive Lagrangian radii closer together without actually indicating a reduction 
in mass-segregation. It is worth noting that mass segregation indeed does become more 
pronounced in DC after 5-10 Trh after the most massive stars have evolved away and the 
mass function has become more peaked. SP and Ca both exhibit extensive mass segregation 
based on the Lagrangian radii with SP exhibiting a slightly larger discrepancy between the 
massive and overall radii than Ca. Again Ca has a shorter two-body relaxation timescale 
than SP and we invoke the same idea as in DC to explain this effect. Thus, since both 
SP and Ca exhibit a similar size at each overall Lagrangian radius, we do not expect two- 
generation clusters to be reliably distinguishable from their single-generation counterparts 
by observations of overall size or mass segregation. 

In Figure H] we consider the final density profiles of all simulated clusters. Note that the 
75% Lagrangian radius is at log r/r vir « 1 and the profile is noisy outside this region due to 
low-number statistics. The density profiles are remarkably similar with all clearly exhibiting 
an exponential halo and a flatter core. Ca and Cb are almost indistinguishable from SP 
except at the very center. DC is slightly less dense than the others but has essentially the 
same profile. The only simulation which displays a significant deviation is Cc with a steeper 
core profile than the other simulations. Cc differs from the other two-generation clusters in 
other ways and will be discussed in more detail in § 14.1.11 

In Figure [5] we consider the radial velocity dispersion profile for all simulated clusters. 
Again the 75% Lagrangian radius is at log r/r vir ps 1 and statistics beyond this are point 
unreliable due to low numbers. We note that the cores of Ca and Cb are slightly more 
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distinct than SP and DC and that all profiles are vertically shifted compared to each other. 
These variations are as large between the two-generation clusters themselves as between the 
two-generation and single-generation clusters and they are merely simulation-to-simulation 
variations. The peaks in Ca and Cb are caused by high- velocity single stars and are a product 
of low-number statistics. The largest outlier is again Cc and we now turn to a discussion of 
this simulation. 



4-1.1. Simulation Cc 

We consider core densities as a fucntion of time for all simulations in Figure O Unlike the 
other simulations, Cc exhibits a clear peak in core density at approximately 15 Tr# 11.5 
Gyr), corresponding nicely to the steeper density profile in the core region in Figure HI We 
tentatively propose that Cc has experienced a core-collapse event at 15 Trh- If correct, 
this is an interesting result and tells us that the dynamics of two-generation clusters are 
similar to those of their single-generation counterparts in yet another way. That only one 
simulation would experience core collapse at a specific time is not surprising since the onset 
of core collapse is associated with individual escapers and hard binaries, the specifics of which 
vary between simulations. We caution, however, that with only 12000 particles such a peak 
could be caused by a few massive stars in the core region rather than a full core-collapse 
event. We also note that the kinetic energy of escapers in a core-collapsed cluster should 
be systematically higher than in a non-core-collapsed cluster due to binary scatting from 
the core. This is not the case for Cc (Table [2]). Further simulations with more particles are 
needed to fully pin down the core-collapse behaviour of two-generation clusters. 

4.2. Stellar Population 

Finally, we consider the stellar populations of the clusters. In Figure [7] we present, 
as histograms, the final mass functions for SP, DC and Ca. In the case of SP there has 
been very little evolution in the mass function. The number of 0.7 — 1M stars has been 
slightly enhanced due to the evolution of high-mass stars into white dwarfs but the stellar 
population is still dominated by low-mass main sequence stars (Tabled]). By contrast both 
DC and Ca are very strongly peaked at 0.7 — 1M . This peak is composed mainly of C-0 
white dwarfs with intermediate-mass progenitors and which make up a much greater fraction 
of the population than they do in SP. Indeed DC is dominated by such compact remnants 
(Table [3]). The strong double peak in the mass function of Ca is also to be noted as it 
resembles the combination of the final mass function for SP and DC, as would be expected 
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for a cluster formed from a combination of two such IMFs. These mass functions are all 
quite different and indeed all three types of clusters should be distinguishable by obtaining 
accur ate white dwarf to main sequence star number ratios. Recent observations of NGC 
6752 (IFerraro et al.ll2003l ) show an anomalously high mass-to-light ratio in the core which 
might be indicative of a large population of low-luminosity, compact objects such as old 
white dwarfs. A early, massive generation could be one way of producing such an excess 
although much more study is needed before we would seriously propose this scenario as an 
option in particular cases. 



5. Conclusions 

We have considered the plausibility of a two-generation formation scenario for galactic 
globular clusters from the standpoint of dynamics and to a lesser extent stellar populations. 
Based on our simulations we expect two-generation clusters to be long-lived objects that 
at the current age of the galactic globular cluster population are very dynamically similar 
to their single-generation counterparts. We tentatively propose that two-generation clusters 
can undergo core collapse at a similar ime to the single generation clusters, although more 
investigation is needed on this front. We expect two-generation clusters with top-heavy IMF 
first-generations to have different stellar populations and specifically we predict that two- 
generation clusters with a top-heavy IMF would be strongly enhanced in white dwarf stars. 
Clusters observed to have such an enhancement would be candidates for the D'Antona and 
Caloi two-generation formation scenario. 

There are two areas in particular that could use further expansion. The first is the 
treatment of the ejecta in the initial generation. In particular a more complete treatment of 
the gas and a full exploration of different star formation efficiencies woul d be useful in order to 



make perdictions about young clusters. In particular it has been shown (IBastian fc Goodwin 



20061 ) that rapid loss of gas from a young cluster can leave the remaining stellar component 
out of equilibrium. The loss of the non-star-forming fraction of the ejecta from the first 
generation could have a similar result on the initial state of a two generation cluster and 
deserves exploration. The second area is to perform further simulations to pin down exactly 
when a two-generation cluster will experience core collapse. 

Based on our simple simulations, we find that two-generation clusters are plausible on 
the basis of d ynamics and co uld form a fraction of the population of galactic globular clusters. 



Interestingly iKarakasi (120061 ) show that AGB stars may not be capable of producing sufficient 
helium enhancement to produce the BHB of NGC 2808 without violating other constraints 
(most importantly C + N + O ~ const.). Thus, although we find two-generation clusters to 
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be a dynamical possibility, the application of the two-generation scenario to the motivating 
example cluster remains in doubt. 



A. Appendix A 



The top-heavy IMF from lDAntona fc Caloil (120041 ) is not included in the standard set 



of Starlab IMF prescriptions and we have encoded it ourselves. First we fix the constants c\ 
and C2 by the number of RHB stars in NGC 2808, the specific choice of Mb, a and j3 and 
the criterion ciM B ^ 1+a ^ = c%M B . This yields: 



dN 
Cl ~dM 



x M~sr ] c 2 = c lM t a) 



M R HB 

We can then calculate the number of stars in the a and /3 regimes: 

—a —p 
The mass for an individual node is then generated using the prescription: 



where P mr and P m are random numbers between zero and one which determine respectively 
the mass regime a star will be in (a or (3) and what its mass in that regime will be. This 
mass-function is defined entirely by the choice of a, (3 and M B and produces the needed 
enhancement in high-mass stars. 
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Fig. 1. — A D'Antona and Caloi IMF (dotted line) prepared by our new Starlab algorithm 
compared with a Salpeter IMF prepared using one of the standard Starlab algorithms (solid 
line) . 
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Fig. 2. — On the left are shown the Lagrangian radii for a cluster with a top-heavy IMF 
evolving for 200 Myr. On the right are the Lagrangian radii for the same cluster evolving in 
a Plummer potential with a static mass equal to the mass lost to stellar evolution over the 
200 Myr life of the cluster. From bottom are the 10%, 50% and 75% Lagrangian radii. Solid 
indicates the entire system, dotted indicates the most massive 10% of stars. 
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Fig. 3. — The Lagrangian radii for a Salpeter IMF (top left), a D'Antona and Caloi IMF 
(top right) and a D'Antona and Caloi IMF with a second generation added according to 
a Salpeter IMF (bottom). From bottom to top in each graph are the 10%, 50% and 75% 
Lagrangian radii. Solid is for all stars, dotted for the most massive 10% of stars. The spike 
just before 10 Trh in SP is created by the ejection of a massive binary from the cluster core. 
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Fig. 4. — The radial density profiles of the simulated clusters. Thin solid is SP, thick solid 
is DC, dotted is Ca, broken is Cb and dot-dashed is Cc. Profiles are taken to the stripping 
radius (100 initial R V i r ). 
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Fig. 5. — Profiles of the radial velocity dispersion for the simulated clusters. Thin solid is 
SP, thick solid is DC, dotted is Ca, broken is Cb and dot-dashed is Cc. Profiles are taken to 
the stripping radius (100 initial R V i r ). 
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Fig. 6. — Core densities for the simulated clusters as a function of time. Thin solid is SP, 
thick solid is DC, dotted is Ca, broken is Cb and dot-dashed is Cc. Time is physical so that 
all simulations can appear on the same graph. 



-22 - 




Fig. 7. — Histograms displaying the final mass function for the clusters SP (top left), DC 
(top right) and Ca (bottom). 
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Simulation Parameters 





1 


st g en 


2 nd gen 


Run ID 


N 


M to t{M Q ) 


N 


M tot {M Q ) 


SP 


12000 


4275 






DC 


12000 


26312 






Ca 


4096 


9093 


8350 


3132 


Cb 


4096 


9053 


8400 


3061 


Cc 


4096 


9032 


8300 


2982 



Table 1: The parameters used to initialize our eight simulations. ID identifies the run. The 
total particle number and total mass are then given for the first and second generation 
respectively (the mass of the second generation matches the mass lost in the first). 



Final Physical Parameters 



ID 


Nf 


M f (M Q ) 


IW(Myr) 


Final Agc(Gyr) 


M lo%m (M ) 


M 10%fm (M Q ) 


K esc (N-body) 


SP 


11189 


2684 


930.2 


12.1 


0.5407 


0.5147 


7.376E-2 


DC 


11507 


7308 


413.9 


11.3 


4.537 


0.8864 


2.307E-4 


Ca 


11718 


4336 


799.8 


12.2 


2.053 


0.7689 


7.650E-3 


Cb 


11842 


4370 


736.3 


12.4 


2.072 


0.7693 


1.764E-3 


Cc 


11661 


4359 


738.3 


12.4 


2.120 


0.7744 


1.831E-3 



Table 2: Column (1) gives the run ID, (2) is the final particle number, (3) is the final mass 
of the system, (4) is the initial half-mass relaxation time, (5) is the final age of the cluster, 
(6) and (7) respectively the initial and final cut-off masses for stars to be included in the 
calculation of the massive Lagrangian radii and (8) is the kinetic energy of all stars with a 
square of velocity greater than four times the RMS velocity. 



-24- 



Final Stellar Populations 



Type 


SP 


DC 


Ca 


Cb 


Cc 


Main Sequence 


10688 


3886 


8736 


8844 


8587 


White Dwarfs 


461 


7453 


2896 


2903 


2997 


Neutron Stars 





30 


15 


16 


11 


Black Holes 


2 


4 


3 


2 


1 


Other 


38 


134 


68 


77 


65 



Table 3: The final stellar populations for each simulation. "Other" includes subgiant, 
Hertzprung gap, horizontal branch etc. 



